Phonons and Conduction in Molecular Quantum Dots: Density Functional 
Calculation of Franck-Condon Emission Rates for bi-fullerenes in External Fields 



o 
o 

O 

Q 



Connie Te-ching Chang, James P. Sethna, Abhay N. Pasupathy, J. Park, D.C. Ralph, P.L. McEuen 
Laboratory of Atomic and Solid State Physics (LASSP), 
Clark Hall, Cornell University, Ithaca, NY 14853-2501, USA 
(Dated: February 6, 2008) 

We report the calculation of various phonon overlaps and their corresponding phonon emission 
probabilities for the problem of an electron tunneling onto and off of the buckyball-dimer molecular 
quantum dot C72, both with and without the influence of an external field. We show that the stretch 
mode of the two balls of the dumbbell couples most strongly to the electronic transition, and in 
turn that a field in the direction of the bond between the two C36 balls is most effective at further 
increasing the phonon emission into the stretch mode. As the field is increased, phonon emission 
increases in probability with an accompanying decrease in probability of the dot remaining in the 
ground vibrational state. We also present a simple model to gauge the effect of molecular size on 
the phonon emission of molecules similar to our C72 molecule, including the experimentally tested 
C140. In our approach we do not assume that the hessians of the molecule are identical for different 
charge states. Our treatment is hence a generalization of the traditional phonon overlap calculations 
for coupled electron-photon transition in solids. 
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I. INTRODUCTION 

Physics is full of examples of phonon-coupled quan- 
tum tunneling events. A classic example from the 1960s 
is the work done with trapped-electron color centers in 
the lattices of the alkali-halidcs 1 . More modern exam- 
ples include the study of how the mobility of intersti- 
tials in metals is modulated by coupling of the defect to 
the resulting distortion of the surrounding lattice^ and 
the study of how the inter-chain hopping by polarons is 
affected by phonon interactions 3 . In these studies and 
others, the frequencies before and after the transition 
were assumed to be unchanged and only the coordinate 
about which the harmonic potential is centered shifts. 
Here, our use of the word phonon, traditionally used for 
plane-wave-like solutions in periodic crystals, for vibra- 
tional normal mode is in the same spirit in this context 
for which we use the term quantum dot, a macroscale 
object, for molecule. 

Over the past several years, several experiments and 
theoretical studies^ 6 - have been done where single 
molecules have been used as the medium for vibration- 
assisted tunneling. Some recent experimental ex- 
amples are measurements done with scanning tunnel 
microscopes^, studies of single hydrogen molecules in 
mechanical break junctions^, and the investigations that 
have directly motivated this work, the three-terminal sin- 
gle molecule transistor experiment a 10 ' 11 where a single 
molecule is deposited between two leads and is subjected 
to both a source-drain and gate bias. This is done in the 
Coulomb blockage regime, where the bias is tuned so that 
sequential transport can occur and a differential conduc- 
tance graph can be plotted. In many of these differential 
conductance graphs, in addition to the main lines due to 
the change in the charge state of the molecule, there are a 
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FIG. 1: C72 with 19 meV stretch mode indicated 



series of sidebands thought to be caused by the coupling 
of the electron to the vibrational modes of the molecule. 

In this paper, we present a general theory for these vi- 
brational overlaps where the vibrational modes of both 
the initial and final electronic states of the molecule are 
considered. Charge dependent hessians and anharmonic 
potentials in the context of single molecule transistors 
have been considered previousl y 12 ' 13 where the molecule 
is assumed to have one dominant mode in each electronic 
state. In the field of chemical spectroscopy, this topic has 
been addressed^ through a general consideration of the 
Franck Condon factors with Duschinsky rotation^ and 
its refinement o 16 ' 17 , which allows for different frequencies 
and eigenvectors between different charged states. Our 
paper considers a realistic model of an N-atom molecule 
with 3N possible modes (for example, the bi-fullerene C72 
with 216 possible modes, see figure ((T|)) and allows the 
calculation of experimental scenarios by combining our 
formulation with results from existing quantum chem- 
istry packages. 

Spectroscopy has long been utilized as a tool in both 



2 



chemistry and physics to study the properties and struc- 
ture of atoms and molecules. Different types of spec- 
troscopy are used for different aims; optical spectroscopy 
for example studies the interaction of electromagnetic ra- 
diation with the sample while this paper addresses dif- 
ferential tunneling spectroscopy. Franck-Condon factors 
serve as a very good tool for analyzing the absorption 
and emission band intensities corresponding to vibra- 
tional levels in atoms and molecules^. Over the years, 
many such molecular vibrational spectra have been cal- 
culated and cataloged using Franck-Condon factor s 19 ' 20 . 
Single molecule transistors offer an opportunity to apply 
the Franck-Condon principles to a new system. Because 
we are dealing with single molecules, we can calculate 
(using quantum chemistry packages) the full vibrational 
profile of both the initial and final electronic states of the 
molecule, and thus calculate the Franck-Condon intensi- 
ties generally. 

Franck was the first to postulate that an electronic 
transition can be accompanied by a vibrational excita- 
tion using classical arguments^. Condon later dupli- 
cated the argument using quantum mechanics and the 
Born-Oppcnhcimer approximation 2 ^. According to Con- 
don, the intensity of a particular transition can be de- 
termined by calculating the transition dipolc moment. 
The electronic component can be factored out, leaving 
the square of the phonon overlap to modulate the total 
transition. 

In section II, we set up our Hamiltonian. IN section III, 
we outline our version of the calculations of the phonon 
overlap calculations iincluding the 'Duschinsky rotation'. 
In section IV, we outline our DFT numerical methods. 
Section VI calculates the zero-field overlaps. Section VII 
addresses the overlaps in a field. Section VIII introduces 
a simple two-ball and spring model for the behavior of the 
stretch mode overlap in the presence of a field. Section 
IX makes contact with dl /dV measurements of the entire 
spectrum and section X concludes. 



Hs = £ T{{c^c d + eld) + £ T^c d + 4cD (4) 

k k 

Here we do not incorportate explicit terms for spin and 
charging effects because we're focused on sequential tun- 
neling events where only one elextron is on the dot at 
a time. Terms like the Coulomb effect in equation © 
below, we disregard in our expression. Equation ^ de- 
scribes the electronic component involving the left and 
right leads and the dot, equation ([4]) describes the tunnel- 
ing component, and equation §5§ are the phonon states. 
Here p is the 3iV-component vector for the momentum 
of the N-atom molecule, M is the mass matrix (diago- 
nal entries giving the masses of the different nuclei in 
groups of three), x is the 37V-dimensional vector for the 
atomic coordinates of the molecule, ri and r 2 are the 
3-/V-dimensional vectors for the minimum energy config- 
urations of the inital and final electronic states, and Ki 
and K2 are the quadratic forms giving the energy near ri 
and r 2 . The 1 and 2 indices reference the charge state 
of the molecule. In our transition, the 1 index refers to 
the molecular state with smaller charge, and the 2 index 
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refers to the higher charge state. T fc is the tunneling 
matrix where the superscripts I and r specifies the left 
or right lead. Finally, c) and c are creation and destruc- 
tion operators, respectively, for electrons. An external 
force f on the atomic coordinates shifts the ground state 
configuration: e.g., r 2 — r 2 + K2 — 1 f. 

Although the phonon states can also be expressed in 
second quantized form via the creation and annihila- 
tion operators for bosonic particles and a, we chose 
to express them in first quantized form to facilitate the 
calculation of the 3N-dimensional overlap integrals be- 
tween different vibrational states of our initial and final 
molecule. 



III. PHONON OVERLAP INTEGRALS 



II. THE HAMILTONIAN 

The Hamiltonian for the molecular dot looks like, in a 
mixed first and second quantized formulation: 



H — Hi + H 2 + H3 
where Hi, H 2 , and H 3 are given by the following: 



#1 



E 



efcC/M- + E £kC k c k + £ dc d c d 



Hi 



+ (1-4*,) [l^x-rO+Kxtx-rO] 
+ 4c4±/2(x-r 2 )tK 2 (x-r 2 )] 



(1) 



(2) 



(3) 



In regarding our problem, we are considering the case 
of incoherent tunneling, where one electron is tunneling 
on or off the molecular quantum dot at a time. We also 
assume that the molecule has time to relax to its min- 
imum energy configuration in between tunneling events 
such that all phonon excitations decay before the next 
tunneling event. In recent experiments 2 ^, long phonon 
lifetimes extending at least fifty times beyond the life- 
times observed in Raman spectroscopy have been mea- 
sured for experiments on suspended carbon nanotubes. 
However, the authors note that the lack of coupling 
to a substrate may account for this increase. In other 
experiments^, the experimental setup was arranged to 
increase the lifetime of the electron compared to the 
phonon, allowing for observation of transient vibronic 
levels. For our calculation, we presume that low currents 
and strong phonon coupling between molecule and leads 
ensures vibrational relaxation between electron tunneling 
events. 
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Secondly, we assume the Born-Oppenheimer approx- 
imation where the total wave function is described by 
\9(z,x)) = \ip x (z)4>(x)} , where x labels the nuclear coor- 
dinates as above, and z labels the electron coordinates. 
Strictly speaking, the ground-state electron wave func- 
tion depends parametrically on the nuclear positions x. 
Based on geometric optimization calculations, we know 
the atomic fluctuations Sx are on the order of picometers, 
so we can safely assume (j> x {z) ~ <j){z) and hence factor 
our wavefunction into a purely electronic component and 
a purely nuclear component. The transition then is given 
by the following matrix elements where T is given by ([2]): 



r /4 - (tf/foaOlTltfifoa;)) 

= ( V l {zW(x)\T\ V f(z)^{x)) 
= {^{z)\T\^(z)){4>f(xW(x)) 



(5) 



Following the Landauer approac h 24 ' 25 i 26 and using 
Fermi's Golden Rule to give us a transition probability, 
we can find the conductance. Sequential tunneling in 
the Coulomb Blockade regime is assumed — for the case 
where the quantum dot has energy levels given by E p 
(For us, E p will represent one of many energy eigenstates 
with both a change in electron charge and the excitation 
of one or more vibrational modes). The formula for the 
conductance is then found to be: 



The linear conductance given in equation ((7|) is de- 
pendent only on the ground state since the excitation of 
phonons is related to the applied bias difference across 
the molecule. As each threshold step in bias is crossed, 
new possible pathways are accessible and the squares of 
their overlaps must be added to the expression. 

This phonon overlap integral is the quantity of inter- 
est since it modulates the total transition rate. Its value 
is a measure of the probability of occurrence of a par- 
ticular transition between the initial vibrational state of 
the initial charge state (assumed to always be the ground 
state) and the final vibrational state of the final charge 
state. This quantity will suppress the total transition 
rate matrix element, leading to less intensity in the line. 
Summing over final states yields l 28 : 



]T|(^(*)|^(x))| 2 = i 



(8) 



with the individual terms representing the probability 
decomposition of the initial state in the eigenstates of the 
final potential. Hence the weight of the original transition 
is spread among the phonon excitations. 



A. Normal modes and phonon wavefunctions 



conductance = 



kT 



p p 



P eq (N)F eq {E p \N) 

p N " P ' " P 

x [1 - f(E p + U(N) - U(N - 1) - E F m 



1 1 T 

where T v is the tunneling rate from the dot electron 
energy level p to the left/right leads, N is the num- 
ber of electrons in the dot before the tunneling event, 
U(N) is the charging energy for N electrons on the dot, 
F eq (E p \N) is the conditional probability in equilibrium 
that the level p is occupied given that the quantum dot 
contains N electrons, Ep is the chemical potential of the 
leads, and P eq (N) is the probability that the quantum 
dot has N electrons in equilibrium. In addition, the elec- 
trons in the leads are in the Fermi distribution f(E). 

Since tunneling rates depend exponentially on separa- 
tion, the tunneling rate through one of the leads, say the 
right one, is often much smaller than the other, with the 

pi pr 

smaller rate acting as the bottleneck, „, p , * ~ T 1 !. Hence 

the dependence of Tp on the final state p determines the 
variation of the conductance with energy. This tunnel- 
ing rate is proportional to |Tfj| , where \Tfi\ is given by 
([5]). For a given electronic transition ip" 1 to <pf , assum- 
ing (at the low temperatures of these experiments) that 
the initial molecular state has no phonon excitations, the 
dependence of this matrix element on the final phonon 
state is given primarily by the phonon overlap integrals 
in equation ([5]), and hence 



G ex Tp oc 



\{<f> f {x)W{x))\ 



(7) 



In 3N dimensions, the phonon Hamiltonian for the ini- 
tial charge state is: 

i?=iptM- 1 p+i(x-r 1 )tK 1 (x-r 1 ) (9) 

For molecules with atoms of unequal mass, transform- 
ing from position space to normal modes becomes much 
simpler if we use the standard trick of rescaling the co- 
ordinates by the square root of the mass and shift the 
origin to ri, the equilibrium configuration of the initial 
charge state: 



Hence: 



y=M 1 /2( x _ ri ). 

ntn i 



Hi = 



(10) 



(11) 



where n = M^P and ft 2 = M _1 / 2 KiM -1 / 2 is a matrix 
with dimensions of frequency squared. 

Similarly, the phonon Hamiltonian for the final charge 
state is: 



Ho = 



)tivr 



i 

lT^ + 2 

1 

2 



_(x-r 2 )tK 2 (x-r 2 ) 



ntn + i(y-A)tfi 2 (y-A) (12) 



where: 



.1/2 



( r 2 - ri) 



(13) 



e F <eV 



is the rescaled atomic displacement due to the change in 
charge state. 
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B. 3N-dimensional wavefunctions and overlaps 

In this section we calculate the transition rate from 
the neutral molecule's ground state to the ground state 
and the various singly excited vibrational states of the 
charged molecule. Our calculation of the Franck-Condon 
factors is thus the one-phonon emission special case of 
the more general Duschinsky rotation calculations in the 
chemistry literature^. We present it here partly be- 
cause we find this special case physically illuminating, 
and partly to introduce our notation. We present in the 
Appendix the more complex calculation of the Franck- 
Condon factor from the neutral ground state to a doubly- 
excited vibrational charged state. 

Using the Hermite polynomials associated with so- 
lutions to the harmonic oscillator (H 1 (x) — 2x and 
H 2 (x) = — 2 + 4a; 2 ), and the expression for the excited 
wavefunctions, we have the 3N-dimcnsional vibrational 
eigenfunctions: 



= iV 2 



N 2 ^Mjh((y- A) -el 2 )) 
xexp(-^(y- A)tft 2 (y- A] 



2h 



N 2 
2^2 



^f 2 (^^W^((y-A). e l 2 ))) 



xexp(- — (y- A)t0 2 (y- A)).(14) 



Here, the A encapsulates the geometric reconfigura- 
tion of the molecule (equation IT5|) . the superscript de- 
notes the initial (1) and final (2) charge states, the first 
subscript is the number of phonons emitted, and the sec- 
ond subscript (if any) is the phonon mode a in which 
they were emitted. The frequency of the phonon mode 
a is given by co a and e« is the orthonormal eigenvector 
of mode a for the molecule in the electronic state i. 

The overlap between the two ground vibrational states 

is 



Oo,o = J dy^l*{y)^l{y) 



dyl NxN 2 exp 1 -y 1 — y 



2h' 



eX p(-(y-A)tg(y-A))}. (15) 



We now rewrite expression (|T5|) so that it contains a 
single gaussian rather than a product of two gaussians: 

N X N 2 / dy[e-^^e-wM-^^(y-*)] 

We want to express the single gaussian as one that is 
centered on a new origin y m ax = A with a new Hessian 



f2 so that the integral is of the form: 

mN 2 J d( y - A) e -^y- A ) 1 ' s (y- A )+ B (17) 

which wc know how to solve. Here B is one of the un- 
knowns we're solving for. It's a constant which will be 
pulled out of the integral with a value given in equation 

3 : 




1 yi 2 

FIG. 2: Wavefunctions 'Pq 1 ' and vPq 2 ', f° r harmonic potentials 
rif and £l| in terms of a two-dimensional rescaled coordinate 
(yi,y2), separated by the rescaled length A = v / m(r2 — r-i). 




yi 

FIG. 3: Overlap integrand corresponding to wavefunctions in 
fig ©, centered on A with quadratic form Cl 2 . 



Setting like quantities equal between expressions ([TBI 
and (117|) . we obtain 

A = (fi 1 +Q 2 )~ 1 n 2 A 
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Q = -(£2i + fi 2 ) 

b = -^(A T n 2 A) + ^A T n 2 (n 1 + n 2 )- 1 n 2 A.(m) 



xe ™- 



A+fi 2 A iA + QA 



Our overlap integral now looks like: 

A^2 / dyle-^y-^'^-^} 



xe 2h 



i r At0 2 Ae57i Atn2<ni + n2) lfi2A 



(19) 



Rewriting the constant part of the integral in terms of 
A and f2, we have: 

N X N 2 J dy [ e -^(y-A) t n(y-A) ]e -^Ata 2 A e iAtQA i 

(20) 

Changing variables to y = y — A, this last integral 
is another multidimensional Gaussian, equaling 1/N 2 , 



where N — y det(^). The ground state to ground state 
overlap is then: 



O n 



N t N: 
N 2 



exp ^A^QA^ exp 



-a.|a). 



The probability of being left in the phonon ground 
state, the tunneling rate T, and the conductance G are 
all suppressed by a factor exp(— g) = |Oq,o| 2 , where 



G = -H\0 ,o\ 2 ). 



(22) 



This defines the total g-factor which we will use to char- 
acterize the overall strength of the phonon coupling. 

We can similarly calculate the overlap between the 
ground initial state and a final state with one phonon 
excited into mode a: 

Oo,i« = / dy^ 1} *(y)*S(y-A) 

= J dy{N ie -^y tn ^N 2 ^JJh((y-A)-e^) 

xex p(-^(y- A ) tfi 2(y- A)) }. 

Combining the exponentials, rewriting them in terms of 
Q and A, we find: 

Oo.ia = N X N 2 J dy{VWS((y-A)-6< 2) ) 

xe -i(y-A)tQ (y -A) e - 5 LAtn 2 A e iA+nAj (2 3 ) 

Changing variables to y = y — A, we have 
O , la = N^ 2 J dy^2uj a /h((y - (A - A)) • e {2) ) 



xe « 



\ 'At 

>e 2fi 



Q 2 A e iAtnA 



N!N 2 J dy{ V^/S(y • ei 2) ) e"^^ 



Since the first term in the last integral is odd in y, it 
must vanish. 

Hence, from equation 1211 the overlap between the 
ground initial state and the excited final state is 



On 



We define: 



O ,o ( y/2w a /H 



\Oo. 



,(2) 



(A -A) 



(24) 



|Oo.o| 2 



AI n 



P, 



ground 



AI, 



ground 



(25) 



which experimentally gives the ratio of the current flow- 
ing emitting one phonon in mode a per electron to the 
current emitting zero phonons (the ratio of the step 
heights in the dl/dV curves). Thus, 



g a = ( VW^ 2) • (A - A)) 



(26) 



In the special case Qi = Q2, where the change in charge 
state does not alter the spring constant matrices Ki and 
K2, the phonon frequencies and normal modes remain 
unchanged. It is well known that the total overlap in- 
tegral is related to the one-phonon emission rates in a 
simple way: specifically G — J2 a 9 a - This is no longer 
the case when the two charge states have different spring 
constant matrices: we must calculate them explicitly^. 
The probability of multiple phonons being emitted into 
distinct phonon modes is given by g a gp ■ ■ ■ |Oo,o| 2 > a s it 
is for the traditionally studied case Qi = Q. 2 . But the 
probability for n phonons to be emitted into the same fi- 
nal state is no longer fj|Oo.o| 2 - We do the calculation of 
two phonons in the Appendix; more general Duschinsky 
rotation calculations can be found in the literature^. 



IV. METHODS 

We used Gaussian 2003, a quantum chemistry package 
to calculate all of the quantities needed in our calculation. 
These quantities include the the force constant matrix K 
for different charge states of the molecule, with dimension 
3iV x 3N. This matrix is related to the f2 2 matrix by the 
equation K — Mil 2 since in the cases of both C140 and 
C72, M commutes with ft. We obtain the vibration fre- 
quency eigenvalues and normal mode eigenvectors from 
fl 2 . 

The program also gives the geometrically minimized 
structures of the molecule for its different charge states r 
and the forces on the atoms / under the influence of an 
external electric field. 

All quantities are calculated under the hybrid B3LYP 
level of theory of the DFT (density functional theory) 
model. The basis set used was STO-3G. The energy 
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and minimum geometric structure were also calculated 
for neutral and charged C72 using the more complete 3- 
21G* basis set. Preliminary comparisons with the more 
complete basis set suggest that qualitative features are 
similar to the simpler basis set. All matrix calculations 
are done under Matlab or its freeware clone GNU Octave. 

Because we were working with molecules of consider- 
able size and were calculating vibrational modes which 
require many electronic relaxation calculations, we used 
the minimal STO-3G basis set for our larger molecule 
(N = 140) and the slightly larger 3-21G* basis set for 
our smaller molecule (N = 72). More complete basis sets 
would capture the polarization and charging effects more 
accurately which would serve to increase our g factors 
since the variation between neutral and charged species 
would be more pronounced. However, our analytic ap- 
proaches and their aim are independent of the details of 
the quantum chemistry calculation. 



V. THE MOLECULES AND THEIR MODES 

Our studies were inspired by work done in the McEuen 
and Ralph groups at Cornell and Berkele y 10 ' 11 ' 29 . Specif- 
ically, we looked at the single molecule transistor made 
up of C14CT 9 -, a molecule whose vibrational modes have 
been modeled and studied experimentally by Raman 
spectroscopy^ -. Ci 40 is comprised of two C70 fullerene 
cages covalently bonded to each other via two C — C 
bonds The dominant mode is the low energy inter-cage 
vibration stretch mode at 11 meV shown schematically 
in figure ((4|. The second molecule studied was based 



TABLE I: Change in distances (Ar) between centers of mass 
of the fullerene cages for C140 during different charge transi- 
tions (Qi — > Q2) where Qi is the initial charge state of the 
molecule and Q2 is the final charge state of the molecule. 
Shown are the results of our DFT simulations, and those of 
our simple model (section |VIII|I . 



Transition Qi — > Q2 


DFT Ar [pm] 


Simple Ar = x[Q2] — x[Qi] 


-v 1 


1.005 


3.16 


1 -» 2 


1.794 


9.26 


2^3 


2.333 


14.8 


3^4 


3.056 


19.5 


4^5 


3.7337 


23.4 



Figures (j4|) and §T§ were produced using Gaussian2003 
to minimize the geometry of the molecule and RASMOL 
to plot out the atom positions. 

For both molecules, the low energy modes correspond 
to large scale motion of the molecules such as the bend- 
ing, twisting, or stretching of the two cages with respect 
to each other (acoustic type vibrations) while higher en- 
ergy modes correspond to motion of the atoms on a 
smaller scale (optical type vibrations). For example the 
15mcV mode corresponds to a see-saw motion of the two 
cages with respect to each other, the 17meV mode corre- 
sponds to a twisting motion of the two cages away from a 
central point, while the higher energy 78 meV mode cor- 
responds to simultaneous deformation of the cages them- 
selves. The vibrations Gaussian calculate are within 5% 
of the experimental values. 




11 meV 



FIG. 4: C140 with stretch mode shown schematically 

on our interest in Cuo- We wanted a molecule with sim- 
ilar properties as C140, but with fewer atoms (C72). The 
aim was to increase the accuracy of the basis set used for 
calculations which would be computationally costly with 
larger molecules. 

Like C140 , the dominant mode indicated in our calcula- 
tions was the inter-cage stretch mode which has an energy 
of 19 meV in C72. The molecule is depicted in figure fTJ). 



VI. BASIC QUANTITIES 

The shift in the geometrically minimized structure of 
the C140 molecule as it acquires an extra electron (charge) 
is the predominant factor in determining the amount of 
phonon emission. If the structure changes little, the over- 
lap between the two ground vibrational states of the ini- 
tial and final charge state of the molecule will be larger, 
which suppresses phonon emission since the overlap is 
a mathematical statement of how likely it is for the 
molecule to remain in the ground vibrational state rather 
than transitioning to a excited vibrational state. 

It is not known what the natural charge states of our 
molecule are on a gold substrate, as used in the exper- 
iments we compare to. A single Cqq molecule typically 
has charge — 2e on gold; doubling this, we anticipate that 
the case of interest may involve a transition from perhaps 
four to five extra electrons on our molecule. 

Table U is a chart of the change in the inter-cage dis- 
tance between the two centers of masses of the fullerene 
cages upon adding an electron. As one can see, the 
distance increment increases as the charges increases. 
Therefore, as the charge on the molecule increases, the 
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0.25 r- 
0.2- 



TABLE II: C72 undergoing different transitions. For conve- 
nience we include columns 4 and 5; their result can be deduced 
from the second and third columns 

Table: Probabilities and g-factors for different transitions 



Transition 


G 


Qa — stretch 


|Oo,o| 2 


\Oo,la=stretch\ 


-» f 


0.960 


0.33 


0.38 


0.125 


1 -> 2 


1.18 


0.406 


0.31 


0.126 


2^3 


1.27 


0.455 


0.28 


0.127 


3^4 


1.29 


0.492 


0.27 


0.135 



0.15- 



8 
60 



0.1- 



0.05 



1 L J .. . 1 .1 nil. 



J . L 



50 



100 150 

Mode [meV] 



200 



0.4 



0.3 




1 


1 1 1 




M S 0.2 










0.1 














III! 


L 1 1. . . 


J ■ _ — 1 



50 100 150 200 

mode (meV) 



FIG. 5: g a for the C72 — ► 1 charge-state transition. The 
large peak is the stretch mode a = 10 at 19 meV. Including 
the two-state emission lines would add an additional peak at 
38 meV (twice the stretch mode) and an otherwise roughly 
continuous background (see Fig. (|10[) . 



molecular incremental distortion Ar increases, and con- 
sequently the probability that the molecule will remain 
in the ground vibrational state after an electron has 
hopped on decreases, leading to stronger phonon side- 
bands in the differential conductance graphs. Table |TT] 
gives for each electronic transition of the molecule up to 
a charge state of 5 extra electrons, the total g-factors 
(equation (|2"2f ) in the absence of an applied field, the 
g-factor associated with the first excited state (eqn I2U)) 
where an inter-cage stretch mode phonon is emitted, the 
probability of the molecule remaining in the ground state 
(|Oo,o| 2 ) 7 and the probability that the molecule's final 
vibrational state is the first excited state of the stretch 

mode {\O Aa=stretch\ 2 )- 

Plotting a graph of the g a factor for the electronic 
transition of a neutral molecule to 1- molecule vs. all 216 
modes (as in figure (O) confirms that the stretch mode 
of the molecule dominates the single phonon emission. 



FIG. 6: g a for the C140 — + 1 charge-state transition. The 
large peak is the a = 10 stretch mode at 11 meV. Again, 
we estimate that the only significant two-phonon line is at 22 
meV. 

We also plot the corresponding graph of g a for C140 in 
figure ©• As the charge state increases, the effects and 
phonon sideband strengths will increase. Two phonon 
emission into separate modes is given by the product of 
their respective single mode emission while two phonon 
emission into the same mode is given by equation (|49[) . In 
any case, the approximation we made for equation (|50p 
suggests that one phonon emission line given by gn me v 
dominates. 

Again, it is the stretch mode (whose identity is con- 
firmed by displacing the equilibrium coordinates of the 
molecule by a distortion that is proportional to the eigen- 
mode) that is important. Two-phonon emission may also 
be significant since experimentallySS, there is sometimes 
a second smaller peak at 22meV which may be due to 
two-phonon emission into the same llmeV mode. Two- 
phonon emission, however, yields a small contribution to 
the conductance. For two phonons emitted into the same 
mode, the contribution is given by the product of the sin- 
gle phonon overlaps. For two phonons emitted into differ- 
ent modes, the contribution cannot be simply described 
by such a product and the complete expression obtained 
from integrating the product of the relevant multidimen- 
sional gaussians is needed. Although we can calculate the 
probability of transition to 2-phonon up to n-phonon vi- 
brational final states, we confine ourselves to one-phonon 
emission in our calculations because as will be illustrated 
in figure (|10[) two-phonon emission contributes a continu- 
ous background with the only sizeable contribution from 
two phonon emission into the llmeV mode. 

VII. CONSIDERING EXTERNAL ELECTRIC 
FIELDS 

In reality, the molecule is not in a vacuum but in a 
real environment of leads and substrate. In the experi- 




FIG. 7: |*| 2 of the HOMO level of C72 under an electric field 
of 4x 10 9 V/m along inter-cage bond, showing the polarization 
of the electron density 



ment o 11 ' 29 , there is a range of g-factors for different ex- 
periments involving the same molecule. This implies that 
environmental effects play an important role and moti- 
vates our calculation of g-factors in the presence of ex- 
ternal holds. We account for one feature of this variable 
environment by applying an electric field to the system. 
This external held can come about as a result of image 
charges that are set up across the substrate or across 
the leads when extra electrons are added to the quantum 
molecular dot. 

In the Gaussian2003 program, we can impose an ex- 
ternal field, relax the electronic wavefunction due to the 
induced polarization and measure the force (expressed as 
a 37V vector, in this case 216-vector) on each atom. The 
external field will polarize the charge on the molecule 
as seen in the following representation in figure |J7J| of 
the highest occupied molecular orbital under the influ- 
ence of an external field along the inter-cage bond of 
the molecule (rendered using the freeware Molden) . This 
force will then act to distort the molecule's atomic con- 
figuration via lattice relaxation, leading to an increase 
in pathways available to the electron via vibration as- 
sisted tunneling. The initial and final configurations 
in equations ©,(G£]),and @ are ri = + K^fi and 

V2 = r , 2°' ) + K 2 " 1 f2, allowing us to calculate the g a factors 
and hence the phonon emission rates from equation (|26[) . 
Figure ([5]) shows that g a for the llmeV line for C140 
increases substantially under an external field. 

As can be seen in the plot, the field does increase the 
g-factor from its bare value. At reasonable fields (those 
that we might expect to find in the experimental litera- 
ture) such as the region where the field « 3 x 10 9 V/m, 
g a for the most represented mode (the stretch mode) in- 
creases to about 0.5. This field would correspond to a 
charge placed 7 angstroms away. And for a field corre- 
sponding to a charge placed 6 angstroms away (the clos- 
est plausible distance), g a becomes around 1.0. However, 
in experiments, the g-factor varies from values of much 




Field [V/m] 

FIG. 8: Field dependence of g a for the C140 11 meV stretch 
mode from the DFT calculation. The solid line is the field de- 
pendence for our simple model calculation which is explained 
further in section IVIIII Experimental values (triangles) are 
taken at zero field, but included on the plot in a vertical col- 
umn for visibility 



less than 1 to values as high as 6. In order to reach these 
quantities in our present theory, we would need to impose 
much higher and unphysical fields. 

Another dependence we examined was the g-factor de- 
pendence of the various modes on the angle of a fixed 
electric field. In figure ([9]), the electric field was fixed at 
a value of 4 x 10 9 V/m. The leftmost figure is the llmeV 
mode - the stretch mode. Following it from left to right 
are the 3.7 meV mode magnified by a factor of 20,000; 
the 2.37 meV mode magnified by a factor of 20; the 15 
meV mode magnified by a factor of 5; the 17 meV mode 
magnified by a factor of 200 and finally the 27.6 meV 
mode magnified by a factor of 500. 

The molecule is oriented such that its long axis is 
aligned vertically. From the figure, we see that there is 
no coupling of the stretch mode (left shape) to the field 
when the field is aligned in a direction perpendicular to 
the stretch mode, and that there is maximum coupling 
in the direction parallel to the direction of the stretch 
mode. Note also that g a is non-zero for a =stretch mode 
even in the absence of a field. The symmetries of the 
plots in figure © reflects the symmetry of the modes 
and how they relate to the symmetry of the applied field. 
C140 has Cih symmetry so it can be generated by a rota- 
tion of angle tt around a fixed axis and a symmetry on a 
plane orthogonal to the fixed axis. The normal modes of 
a molecule also possess a definite symmetry with respect 
to the planes of symmetry of the molecule. The symme- 
try of the stretch mode is even under reflection in the x-y 
plane, coinciding with the symmetry of the field E x and 
orthogonal to the field E z . Thus, the strongest coupling 
of the stretch mode is to E x . 
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FIG. 9: Angle Dependence of g a = stretch for C140. Leftmost 
figure is g plotted as function of the angle of electric field 
(an extremely high field magnitude of 5 x 10 12 V/m) for the 
llmeV stretch mode; remaining figures to the right are for 
other modes, and have been magnified considerably to show 
their shape. The vertical represents fields along the long axis 
of the molecule. 



VIII. A SIMPLE MODEL FOR OVERLAPS AND 
FIELDS 

To what extent are these quantum overlaps a result of 
complex quantum chemistry (bonding and anti-bonding 
and electronic rearrangements inside the two cages)? 
How much can we understand from simple electrostat- 
ics of dumbbells? By modeling the system simply as 
two rigid balls connected by a spring subject to an exter- 
nal field, we can obtain some understanding at how the 
dimensions of the problem as well as simple quantities 
might affect the overlap and g- factor. 

We write down the total energy of the system and then 
minimize the energy with respect to the parameters of 
our problem and in the presence of an external field - for 
our case we choose to minimize the charge on one ball 
and the distance x between the two balls. 

The quantities we take into account are as follows: 



E, 



^spring 
Efi e ld 

Ecoulomb 
capacitance 



-K(x 2 - X\ 

q\Ex\ + q 2 Ex 2 

Kqiq 2 
(x 2 - Xi 

lM + _ 

2 r 2 



lkqi 



(27) 



where a is the equilibrium distance of the spring, x\ and 
x 2 are the coordinates of the two balls, r is their radius, K 
is the spring constant of the system, and k is the Coulomb 
constant. 

We also note that M total = M baUl + M baU2 = 2M baU 
and M red = m^+m,'^ = M baU /2 are the well-known 
center of mass and reduced mass for the system. The 
last assignment we make are expressions for the charges 
on each ball {q\ and q 2 ) in terms of the charges in the 
system: 



?2 



9. + 1 

2 2 

Q_q 
2 2 



where Q is the total charge of the system and q is the 
difference between the charges on the two balls. The 
potential energy U then becomes: 

U E S p r i TL g -\- E field ~t~ E/coulomb ~t~ E ca p ac itance 

= -j * , { -2a 2 Eqr + K(q 2 (x-r)+Q 2 (x + r)) 
4r(a + x) I 

+ 2rx(2EQX - Eqx + kx 2 ) 

+ a[K(q 2 + Q 2 ) + 2r(2EQX - 2Eqx + kx 2 )}}. (29) 



Here x = x 2 - 
balls and X 



■ x\ is the relative distance between the two 
= Xl + X2 is the center of mass coordinates 
of the system. We next take the derivative of the poten- 
tial with respect to q the difference in charges on the two 
balls and set the resulting expression (^) equal to zero. 
Solving this expression for q gives us the minimized dis- 
tribution of charges on the balls under an external field: 



Er{a- 



K(a — r + x) 



(30) 



Similarly, we take the derivative of the potential energy 
with respect to the deviation from equilibrium x (^7), 
set this expression equal to zero, and solve for x. We 
keep terms up to second order in Q and E and get: 



x[Q] = AE 2 + BQ 2 + CE 2 Q 2 
where A, B, and C are given by: 



(31) 



.4 
B 
C 



2a 2 - bar - 3r 2 
4kK(a- r) 3 



K(a 2 



2ar) 



Aa 2 k(a - 
2r-a 
8k 2 {a 



r) 3 



(32) 



In table (jT|, we compare the Ar = x[Q 2 ] — x[Qi] in the 
absence of a field for our the simple model and the full 
DFT calculation discussed earlier where Q\ is the total 
charge for the initial system and Q 2 is the total charge 
for the final system. The simple model has between three 
and six times the distortion of the quantum chemistry 
calculation, likely due to a combination of more effective 
screening of the Coulomb repulsion between cages and 
quantum chemistry effects. 

The constants from equation [32] which define the ex- 
pression for x in equation (|31|> . in combination with the 

formula for the zero-point motion xq — \J M f daJ ~ an( l the 

one-dimensional equivalent of equation (|26[) gives us an 
expression for G: 



g = G = 



(xjQ^-xjQ^ 
Ax 2 



(33) 



where in this one-mode limit the log G of the total overlap 
equals the one-phonon emission ratio g. Here M red is is 



equal to 



70mc.rti 



(28) 

point motion xq for C140 is 2.17 pm 



— , ujq = lu stretch = llmeV and the zero- 
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Therefore, using these formulas our g factor for the 
— * 1 transition of Cu is 0.535 (and for C72 it is 0.92). 
The complete 3iV-dimensional calculations in the absence 
of a field for the same transition yields a g a for the stretch 
mode of 0.23 for C 140 (and 0.33 for C* 72 ). This differ- 
ence is not as large as one would expect from the differ- 
ence in the center-of-mass motions: the llmeV stretch 
mode incorporates motions that do not simply change 
the center of mass separation. In total, the many-body 
DFT calculations show a stretch-mode phonon emission 
about a factor of three smaller than that predicted from 
the simple physical model, and captures the cage size- 
dependence of g, rather well. 

Finally, we compare the field dependence of the g a be- 
tween the simple model and the DFT calculation given 
in figure ||5J). The field dependence works out quite well. 
This simple model could be made more realistic by in- 
corporating features from the DFT calculation (such as 
charging energies), but that would take us beyond our 
current illustrative goal. 
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FIG. 10: I-V curve predicted for C72 for one-phonon process 
(solid line) and up to two-phonon processes (approximate, 
dashed line), using the DFT STO-3G basis set. The arrow 
indicates the position of the two-phonon contribution from 
the stretch mode. 



IX. CURRENT DUE TO PHONON 
TRANSITIONS 



tions [34l and [ 



0O\2a 



-GJ2 



9 l M 



(34) 



Using the g-factors corresponding to all of the different 
single phonon modes, we graphed a current vs. volt- 
age graph for C72 using the simplified formula where 
all the phonons are identical in both charge states of 
the molecule. Figure (11) gives the current divided by 
Iq vs. the available energy above the ground state to 
ground state threshold for both one-phonon emission 
processes (solid line) and up to two-phonon processes 
(dashed line). The plots are constructed by iteratively 
calculating phonon emission from a pool of available en- 
ergy. As energy decreases, less is available for emitting 
phonons. Our g a 's make use of the fact that the phonon 
quadratic forms change between different charge states. 
As you can see, the currents due to one-phonon processes 
and for up to two phonon processes share similar gross 
features at the beginning such as the jump in current at 
the 19 meV energy mode corresponding to the stretch 
mode of the molecule. However they start to deviate as 
energy increases until they level off at different values of 
current (0.8 for the one-phonon process and 0.95 for the 
two-phonon processes) which would seem to indicate that 
two-phonon processes will play a role in the I-V charac- 
teristics of a molecular quantum dot. 

In addition, the I-V curve that includes all n-phonon 
processes will asymptote to one. The two-phonon con- 
tribution forms almost a continuous background, except 
for 2u stretch, whose position is shown with an arrow in 
figure (11). We also note that the our treatment of 
two-identical-phonon emission is (for convenience) not 
the correct formula derived in equation (|49[) which al- 
lows the frequencies to change between the initial and 
final states, but the approximate formulas given by equa- 



I/I (E) = ]T(0 

a 

+ ^ O 0jlala >@(E - Kv a - Kv a ') (35) 

As one may observe, the approximate two-phonon rates 
form a featureless background except for the 2-stretch- 
mode phonon peak, a result of the weak couplings to the 
other modes. 

X. CONCLUSION 

There is much recent interest in vibrating mechanical 
systems coupled to electron transport on the nanoscale, 
from nanomechanical resonators 31 ' 32 to single-electron 
shuttle s 33 ' 34 . Vibrational effects on electron transport 
through molecules have been studied since the 1960s in 
devices containing many molecules 3 -^, and more recently 
have been shown to be important in transport through 
single molecules measured using scanning tunneling mi- 
croscopes^, single-molecule transistors 1 ^ 1 -, and mechan- 
ical break junctions^. In a natural extension of work 
done in the 1920s by Franck, Condon, et al. in atomic 
spectra, we have studied the effects of molecular vibra- 
tions on electron transport through a molecule. We have 
shown that density functional theory calculations of the 
normal modes and deformations, coupled to a straightfor- 
ward linear algebra calculation, can provide quantitative 
predictions for the entire differential tunneling spectrum, 
even including external fields from the molecular environ- 
ment. 
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XI. APPENDIX 

Here we show how one can calculate the Franck- 
Condon factor for a transition from the neutral ground 
state to an excited state with one vibrational mode in 
a doubly-excited state. (For emission into general ex- 
cited states, we would need to use the appropriate multi- 
dimensional gaussian multiplied by the appropriate Her- 
mite polynomials. This calculation quickly becomes com- 
plicated^, and for the molecules of interest to us, mul- 
tiple phonon emission is rare.) From the vibrational 
states given in (fl"4"|) , we're interested in the following 
integral: 



o ,2 = J dy*o 1} *(y)^(y- 



(36) 



where we can split the integral into two parts: 

= j dyN.N^^ ■ (y- A)] 2 e^^( yt ° iy ) 
xe -A((y~ A ) tQ 2(y-A)) 
(37) 

(38) 

Expression ([38f is just — ^Oo.o'i we concentrate on 
expression (l37|) . First, as we did for the Oo,i case, we 
rewrite this integral in terms of the quantities A and Q: 

A^ 2 V2(^) J dy[ zg) . (y _ A)] 2 e -*(y-A)tfi (y -A) 

(39) 



xe-* Atn2A e* A ( ni+n2 ) A 



We want to make this expression look like the known 
gaussian integral: C\ J dxx 2 e~ x +C2 where C\ and Ci 
are constants. Changing variables to y: 

y =_y- A 

dy = dy 

y = y + A (40) 
we rewrite the integral as one over d n y: 

Nl N 2 V2(^) f df[eW ■ (y + A - A)] 2 e "^ 



Expanding out the term in brackets in equation (141|) , we 
get: 



e( 2 ).(A-A)] 2 = [ei 2 ).y + d 2 ] 2 
(ei 2 ).y) 2 + 2d(6( 2 ).y) + d 2 (42) 



where d = e^ 2) ■ (A - A). 

The second term in equation (142|) will be zero in the 
integral because of symmetry considerations which dic- 
tate that odd powered gaussian integrals of the form: 
J dxx n e~ x where n is odd always equal zero.. The only 
terms in the integral of equation (|41[) that remain are the 
first term and the constant d 2 . 

We transform this integral into the appropriate normal 
mode basis. Since, we are integrating over the coordi- 
nates centered on A for a system with a Hessian of Q, 
we want to rewrite everything in terms of the eigenmodes 
of the averaged Q. We'll call these eigenmodes pp where 
the following definitions hold: 



y = J2 p PP0 





(43) 



Here, pp are the orthonormal eigenvectors for Q. and 
P/3 are the weightings of each mode's contribution to y. 
Hence: 

Again, the second term is odd in the new integration 
variables pp and will be zero. Rewriting the integral in 
dp and remembering that p diagonalizes f2, the integral 
from equation ([36]) and (|4Tj) becomes: 



xe -5k At " 2A e^ A ( ni+n2 ) A 



P J 
xe^57r Atn2A e^ A(ni+n2)A . 



But fx 2 e- Ax *dx = = ^fe~ At 'dx i so 



(45) 



d n 



2 „4Ey p^v 



1 



h l 



I irft, -p-i- / irft 
up -f-L y u g i 



(46) 
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Hence, the first term in (|42]) from equation (|45|) becomes: 

(47) 

which from equation (|2ip we see is 

Oo,oV2^E^(^ 2) -^) 2 - ( 48 ) 

Combining this with the third term from equation (|4"2"|) 
and expression ([38]) , our expression for the 0-^2 overlap 
becomes: 

+ V2(^)[e^ ■ (A - A)] 2 - ±} 

+ V2(^)[e^-(A-A)r-±}. (49) 



If Q.I = ft 2 (i-e., there is no change in the harmonic 
potential), ft = ft2 and hence e^-pp — 5 a i3 and uip = uip. 
The first sum reduces to -^=, canceling the last term. 
Therefore in this case: 

|^| 2 = 2& 2 [^.(A-A)f = f. (50) 



The change in harmonic potential upon charging the 
molecule allows for phonon emission even in the absence 
of a configurational shift. Therefore, even if A = A = 0, 
phonons can be emitted both because of frequency shifts 
(uj a ^ tip) and because the normal modes change (e Q ^ 
Pp)- 

The other overlaps could be performed in a similar way. 
The strategy is to write everything in terms of integrals 
of J dxx n e~ x by transforming to the basis of the aver- 
aged gaussian with ft. The odd powered integrals are 
eliminated and the even powered terms remain. 
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